//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : WavesComputer.h
//----------------------------------------------------------------------------

#ifndef WAVESCOMPUTER_H
#define WAVESCOMPUTER_H

#include "defines.h"
#include "matvec.h"
#include "RrbAkashSolver.h"

#ifdef USE_CUDAWAVESCOMPUTER
#include "CudaWavesComputer.h"
#endif

class Solver;
class LowPassFilter;
class WeakRefl;
class IncomingWaves;


/// computes waves according to Klopman's model
class WavesComputer
{
public:

    WavesComputer(int Nx1, ///< number of data points in x1 direction
                  int Nx2, ///< number of data points in x2 direction
                  const Array2D<REAL>& h0, ///< array with depth data
                  const Array2D<REAL>& U1, ///< array with speed data in x1 direction
                  const Array2D<REAL>& U2, ///< array with speed data in x2 direction
                  const Array2D<SHORT>& wet, ///< array with wet/dry data
                  Array2D<REAL>* pzeta, ///< array with waveheight data
                  Array2D<REAL>* pzetam1, ///< array with previous waveheight data
                  Array2D<REAL>* pphi, ///< array with speed potential data
                  Array2D<REAL>* pphim1, ///< array with previous speed potential data
                  Initialization init, ///< perform wave initialisation?
                  REAL time, ///< start time if not zero
                  REAL deltat, ///< last timestep if not zero start time
                  REAL pstart, ///<< [s] fade-in time of pulses and incoming waves
                  const char *solvername, ///< name of solver to use
                  REAL dx1, ///< grid size in x1 direction
                  REAL dx2, ///< grid size in x2 direction
                  REAL g, ///< gravity
                  REAL tau_brk, ///< wave breaking parameter
                  int wvin_spec_count, ///< number of specifications for incoming waves
                  const Array1D<REAL> &wvin_spec_ampl, ///< incoming waves amplitudes
                  const Array1D<REAL> &wvin_spec_phase, ///< incoming waves phases
                  const Array1D<REAL> &wvin_spec_freq, ///< incoming waves frequencies
                  const Array1D<REAL> &wvin_spec_dir, ///< incoming waves directions
                  int bcN, ///< boundary condition north
                  int bcE, ///< boundary condition east
                  int bcS, ///< boundary condition south
                  int bcW, ///< boundary condition west
                  int somfeld, ///< absorb. bc type: 1=classical Sommerfeld, 2=improved
                  REAL k_wrbc, ///< weak reflection parameter
                  VariationalModel varmod, ///< variational model type
                  REAL kappa_c, ///< wave number
                  REAL psitol, ///< iterative solver tolerance
                  bool short_wave, ///< short wave equations?
                  bool upwind, ///< upwind scheme?
                  REAL eps, ///< smoothing tolerance
                  int lowpass_halfwidth); ///< low pass filter half width

    virtual ~WavesComputer();

    void compute(REAL dt, ///< delta T
                 const Array2D<REAL> &pressurePulse, ///< array with pulse pattern
                 double *pCgTime); ///< time used by solver

    bool save(const char* filename);
    bool load(const char* ffilename);

    REAL getTime() const;
    REAL getDeltaTime() const;

private:

    /// @cond nevercall
    // should never be called!
    WavesComputer();
    WavesComputer(const WavesComputer& m);
    WavesComputer & operator =(const WavesComputer& m);
    /// @endcond

    void initialize(const char *solvername, ///< name of solver to use
                    int somfeld, ///< absorb. bc type: 1=classical Sommerfeld, 2=improved
                    REAL k_wrbc, //
                    VariationalModel varmod, //
                    REAL kappa_c, //
                    REAL psitol, /// iterative solver tolerance
                    int wvin_spec_count, //
                    const Array1D<REAL> &wvin_spec_ampl, // 
                    const Array1D<REAL> &wvin_spec_phase, //
                    const Array1D<REAL> &wvin_spec_freq, // 
                    const Array1D<REAL> &wvin_spec_dir, //
                    Initialization init); ///< initialize waves?

    void startup();
    void waveBreaking();
    void boundaryZetaPhi();
    void solvePsi(REAL pslow);
    void boundaryPsi();
    void weaklyReflective();
    void leapFrog();
    void timeDerivatives(REAL pslow, const Array2D<REAL> &pressurePulse);
    void coupledMatrices();
    void shiftZetaPhi();
    void shiftPsi();
    void adjustPsi();
    void smoothZetaPhi();
    void constructBoundaryStencil();
    void storeRHS();
    void ComputeRHS(Array2D<REAL> &RHS_ZZ, Array2D<REAL> &RHS_FF, REAL pslow, const Array2D<REAL> &pressurePulse);
    void ComputeRHS(Array2D<REAL> &RHS_PP, REAL pslow);

    int m_Nx1;
    int m_Nx2;
    REAL m_dx1;
    REAL m_dx2;
    REAL m_dx1square;
    REAL m_dx2square;
    REAL m_dx1twice;
    REAL m_dx2twice;
    REAL m_g;
    REAL m_tau_brk;
    int m_bcN;
    int m_bcE;
    int m_bcS;
    int m_bcW;
    bool m_short_wave;
    REAL m_eps;
    bool m_upwind;
    int m_lp;
    const Array2D<REAL> &m_h0;
    const Array2D<REAL> &m_U1;
    const Array2D<REAL> &m_U2;
    const Array2D<SHORT> &m_wet;
    Array2D<REAL>& m_zeta;
    Array2D<REAL>& m_zetam1;
    Array2D<REAL>& m_phi;
    Array2D<REAL>& m_phim1;
    REAL m_pstart;
    REAL m_time;
    REAL m_dt_old;
    REAL m_dt;

    Solver *m_pSolver;
    RrbAkashSolver *m_zfSolver;
    IncomingWaves *m_pIncomingWaves;
    WeakRefl *m_pWrBc;
    LowPassFilter *m_pLpFil;

    REAL m_pi;
    REAL m_dt_new;
    REAL m_zeta_mean;
    REAL m_sum_en_wgt;
    Array2D<REAL> m_zetam2;
    Array2D<REAL> m_phim2;
    Array2D<REAL> m_psi;
    Array2D<REAL> m_psim1;
    Array2D<REAL> m_psim2;
    Array2D<REAL> m_dt_phi;
    Array2D<REAL> m_dt_phim1;
    Array2D<REAL> m_dt_zeta;
    Array2D<REAL> m_dt_zetam1;
    Array2D<REAL> m_matS;
    Array2D<REAL> m_matC;
    Array2D<REAL> m_matW;
/* Members added by Akash*/
    Array2D<REAL> m_matE;
    Array2D<REAL> m_matN;

// Matrix Diagonals for Zeta -- Only three are stored as the matrix is skew symmetric.
    Array2D<REAL> m_matzS;
    Array2D<REAL> m_matzW;
    Array2D<REAL> m_matzC;
    Array2D<REAL> m_matzN;
    Array2D<REAL> m_matzE;


// Matrix Diagonals for phi -- Only three are stored as the matrix is skew symmetric.

    Array2D<REAL> m_matfS;
    Array2D<REAL> m_matfW;
    Array2D<REAL> m_matfC;
    Array2D<REAL> m_matfE;
    Array2D<REAL> m_matfN;


// Matric for coupling
    Array2D<REAL> m_matzfS;
    Array2D<REAL> m_matzfW;
    Array2D<REAL> m_matzfC;
    Array2D<REAL> m_matzfN;
    Array2D<REAL> m_matzfE;

    Array2D<REAL> m_matfzS;
    Array2D<REAL> m_matfzW;
    Array2D<REAL> m_matfzC;
    Array2D<REAL> m_matfzN;
    Array2D<REAL> m_matfzE;


    Array2D<REAL> m_matzpS;
    Array2D<REAL> m_matzpW;
    Array2D<REAL> m_matzpC;
    Array2D<REAL> m_matzpN;
    Array2D<REAL> m_matzpE;

    Array2D<REAL> m_RHSzz;
    Array2D<REAL> m_RHSff;
    Array2D<REAL> m_RHSpp;


    std::vector<Array2D<REAL> > m_matzz;
    std::vector<Array2D<REAL> > m_matff;
    std::vector<Array2D<REAL> > m_matzf;
    std::vector<Array2D<REAL> > m_matfz;
    std::vector<Array2D<REAL> > m_matzp;
    std::vector<Array2D<REAL> > m_matpf;





// --------------------------------------------------------------- //
    Array2D<REAL> m_bpsi;
    Array2D<REAL> m_alfa_brk;
    Array2D<REAL> m_zeta_std;
    Array2D<REAL> m_SffC;
    Array2D<REAL> m_SRCp;
    Array2D<REAL> m_SRCz;
    Array2D<REAL> m_SppC;
    Array2D<REAL> m_SzzC;
    Array2D<REAL> m_SzpS;
    Array2D<REAL> m_vD0;
    Array2D<REAL> m_en_wgt;
    // Boundary Node variables ------------------------ //
    Array1D<REAL> m_cN;
    Array1D<REAL> m_cS;
    Array1D<REAL> m_stoN;
    Array1D<REAL> m_stoS;
    Array1D<REAL> m_stozzN;
    Array1D<REAL> m_stozzS;
    Array1D<REAL> m_stozfN;
    Array1D<REAL> m_stozfS;
    Array1D<REAL> m_stoffN;
    Array1D<REAL> m_stoffS;
    Array1D<REAL> m_zeta_in_N_ext;
    Array1D<REAL> m_zeta_in_S_ext;
    Array1D<REAL> m_zeta_in_N_cen;
    Array1D<REAL> m_zeta_in_S_cen;
    Array1D<REAL> m_zeta_in_N_int;
    Array1D<REAL> m_zeta_in_S_int;
    Array1D<REAL> m_phi_in_N_ext;
    Array1D<REAL> m_phi_in_S_ext;
    Array1D<REAL> m_phi_in_N_cen;
    Array1D<REAL> m_phi_in_S_cen;
    Array1D<REAL> m_phi_in_N_int;
    Array1D<REAL> m_phi_in_S_int;
    Array1D<REAL> m_psi_in_N_ext;
    Array1D<REAL> m_psi_in_S_ext;
    Array1D<REAL> m_psi_in_N_cen;
    Array1D<REAL> m_psi_in_S_cen;
    Array1D<REAL> m_psi_in_N_int;
    Array1D<REAL> m_psi_in_S_int;
    Array1D<REAL> m_cE;
    Array1D<REAL> m_cW;
    Array1D<REAL> m_stoE;
    Array1D<REAL> m_stoW;
    Array1D<REAL> m_stozzE;
    Array1D<REAL> m_stozzW;
    Array1D<REAL> m_stozfE;
    Array1D<REAL> m_stozfW;
    Array1D<REAL> m_stoffE;
    Array1D<REAL> m_stoffW;
    Array1D<REAL> m_zeta_in_E_ext;
    Array1D<REAL> m_zeta_in_W_ext;
    Array1D<REAL> m_zeta_in_E_cen;
    Array1D<REAL> m_zeta_in_W_cen;
    Array1D<REAL> m_zeta_in_E_int;
    Array1D<REAL> m_zeta_in_W_int;
    Array1D<REAL> m_phi_in_E_ext;
    Array1D<REAL> m_phi_in_W_ext;
    Array1D<REAL> m_phi_in_E_cen;
    Array1D<REAL> m_phi_in_W_cen;
    Array1D<REAL> m_phi_in_E_int;
    Array1D<REAL> m_phi_in_W_int;
    Array1D<REAL> m_psi_in_E_ext;
    Array1D<REAL> m_psi_in_W_ext;
    Array1D<REAL> m_psi_in_E_cen;
    Array1D<REAL> m_psi_in_W_cen;
    Array1D<REAL> m_psi_in_E_int;
    Array1D<REAL> m_psi_in_W_int;
    Array2D<REAL> m_phiN;
    Array2D<REAL> m_psiN;
    Array2D<REAL> m_phiS;
    Array2D<REAL> m_psiS;
    Array2D<REAL> m_QphiN;
    Array2D<REAL> m_QphiS;
    Array2D<REAL> m_QpsiN;
    Array2D<REAL> m_QpsiS;
    Array2D<REAL> m_QzetaN;
    Array2D<REAL> m_QzetaS;
    Array2D<REAL> m_zetaN;
    Array2D<REAL> m_zetaS;
    Array2D<REAL> m_phiE;
    Array2D<REAL> m_psiE;
    Array2D<REAL> m_phiW;
    Array2D<REAL> m_psiW;
    Array2D<REAL> m_QphiE;
    Array2D<REAL> m_QpsiE;
    Array2D<REAL> m_QphiW;
    Array2D<REAL> m_QpsiW;
    Array2D<REAL> m_QzetaE;
    Array2D<REAL> m_QzetaW;
    Array2D<REAL> m_zetaE;
    Array2D<REAL> m_zetaW;
    // -------------------------------------- //

#ifdef USE_CUDAWAVESCOMPUTER
    CudaWavesComputer *m_pCudaWavesComputer;
#endif
    
};

#endif // WAVESCOMPUTER_H
